Hydrogen and muonium in diamond: A path-integral molecular dynamics simulation 
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Isolated hydrogen, deuterium, and muonium in diamond have been studied by path-integral molec- 
ular dynamics simulations in the canonical ensemble. Finite-temperature properties of these point 
defects were analyzed in the range from 100 to 800 K. Interatomic interactions were modeled by a 
tight-binding potential fitted to density-functional calculations. The most stable position for these 
hydrogenic impurities is found at the C-C bond center. Vibrational frequencies have been obtained 
from a linear-response approach, based on correlations of atom displacements at finite temperatures. 
The results show a large anharmonic effect in impurity vibrations at the bond center site, which 
hardens the vibrational modes with respect to a harmonic approximation. Zero-point motion causes 
an appreciable shift of the defect level in the electronic gap, as a consequence of electron-phonon 
interaction. This defect level goes down by 70 meV when replacing hydrogen by muonium. 

PACS numbers: 61.72.-y, 71.55.Cn, 81.05.Uw 
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I. INTRODUCTION 

In spite of being one of the simplest possible impu- 
rities, a detailed understanding of the physical proper- 
ties of isolated hydrogen in group-IV materials requires 
the combination of advanced experimental and theoret- 
ical methodsi^ Experimental investigations on atomic, 
interstitial hydrogen have been so far scarce, since H has 
turned out to be difficult to observe as an isolated im- 
purity in these materials. In silicon, electron paramag- 
netic resonance experiments^ indicated that hydrogen 
is located at a bond-center (BC) site, midway between 
two nearest-neighbor host atoms, in agreement with sev- 
eral theoretical calculationsiSi^ However, for isolated hy- 
drogen in diamond there is little detailed experimental 
information available in the literature. In this context, 
muonium may be considered as a light pseudoisotope of 
hydrogen, due to the small mass of the muon /i+ (about 
1/9 that of the proton). This means that in many re- 
spects muonium is expected to behave similarly to hydro- 
gen, with the appropriate corrections due to zero-point 
motion and related effects. Thus, muon implantation 
experiments in diamond as well as several theoretical ap- 
proaches have shown that this impurity is metastable at 
a tetrahedral interstitial site (the so-called 'normal muo- 
nium'), and has its lowest energy at or around the bond 
center site ('anomalous muonium' )^iSii 

Apart from its basic interest as an isolated impurity, 
an important property of hydrogen in solids is its ability 
to form complexes and passivate defects. This has been 
extensively studied in the last twenty years for silicon and 
germanium;!'^ More recently, it has been noted that hy- 
drogen can also passivate impurities in diamond.*'^ Over 
the last years, there has been an intense theoretical activ- 
ity on hydrogen-related defects in this material«i24ii*iSii^ 



For isolated hydrogen in the neutral charge state, most 
of the calculations found the BC site as the lowest-energy 
position for this impurity in diamond. Some deviations 
from the BC site were found by Saada et al}^ from a 
tight-binding (TB) calculation. Also, for hydrogen in the 
positive charge state, density-functional theory calcula- 
tions have found this impurity to be stable off-axis in a 
buckled bond-centered configuration, l^ii^ 

In standard electronic-structure calculations in con- 
densed matter, despite of their quantum mechanical char- 
acter, atomic nuclei are treated as classical particles, and 
therefore typical quantum effects like zero-point vibra- 
tions are not directly accessible. These quantum ef- 
fects can be important for vibrational and electronic 
properties of light impurities like hydrogen, especially 
at low temperatures. In fact, the importance of taking 
into account the quantum character of proton and muon 
to study hydrogen-related point defects in diamond has 
been emphasized by Kerridge et ali^ In particular, quan- 
tum tunneling is relevant to understand the properties of 
vacancy-hydrogen complexes in this material. 

Finite-temperature properties of hydrogen-related de- 
fects in solids have been studied by ab initio and tight- 
binding molecular dynamics (MD) simulations. In many 
previous applications of these methods, the atomic nu- 
clei were treated as classical particlesii^iSSiSi In order 
to consider the quantum nature of the nuclei, the path- 
integral molecular dynamics (or Monte Carlo) approach 
has proved to be very useful. A remarkable advantage 
of this method is that all nuclear degrees of freedom can 
be quantized in an efficient manner, thus permitting the 
inclusion of both quantum and thermal fluctuations in 
many-body systems at finite temperatures. In this way, 
Monte Carlo or molecular dynamics sampling applied to 
evaluate finite-temperature path integrals allows one to 
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carry out quantitative and nonperturbative studies of 
highly-anharmonic effects in solidsi^^iS^ 

In this paper, the path- integral molecular dynamics 
(PIMD) method is used to investigate the role of the 
impurity mass on the properties of hydrogenic point de- 
fects. We consider isolated hydrogen, deuterium (D) and 
muonium (Mu) in diamond, in their neutral charge state. 
Special attention has been paid to the vibrational proper- 
ties of these impurities, by considering anharmonic effects 
on their quantum dynamics. The results of the present 
calculations show that such anharmonic effects lead to a 
significant deviation of the vibrational frequencies of the 
impurities, as compared to a harmonic approximation. 
Also, zero-point motion causes an appreciable shift of the 
defect levels in the electronic gap, which is a manifesta- 
tion of the electron-phonon interaction induced by the 
hydrogenic impurities. Path-integral methods analogous 
to that employed in this work have been applied earlier to 
study hydrogen in metak^^'^^ and semiconductors iSSiSSiSl 

The paper is organized as follows. In Sec. II, we de- 
scribe the computational method and the models em- 
ployed in our calculations. Our results are presented in 
Sec. Ill, dealing with the energy of the defects, vibra- 
tional frequencies, and defect levels in the electronic gap 
of diamond. Sec. IV includes a discussion of the results 
and a summary. 



II. COMPUTATIONAL METHOD 
A. Path-integral molecular dynamics 

In the path-integral formulation of statistical mechan- 
ics, the partition function is evaluated through a dis- 
cretization of the density matrix along cyclic paths, 
composed of a finite number L (Trotter number) of 
"imaginary-time" steps^SiSS In the numerical simula- 
tions, this discretization gives rise to the appearance of 
L "beads" for each quantum particle. Thus, this method 
exploits the fact that the partition function of a quan- 
tum system is formally equivalent to that of a classical 
one, obtained by replacing each quantum particle by a 
ring polymer consisting of L particles, connected by har- 
monic springs>22i2^ In many-body problems, the configu- 
ration space is usually sampled by Monte Carlo or molec- 
ular dynamics techniques. Here, we have employed the 
PIMD method, which has been found to require less com- 
puter time resources in our problem. Effective algorithms 
to perform PIMD simulations in the canonical NVT en- 
semble have been described in detail by Martyna et alm^ 
and Tuckermani^^ All calculations presented here were 
carried out in the canonical ensemble, using an originally 
developed MD software, which enables efficient PIMD 
simulations on parallel supercomputers. 

The calculations have been performed within the adia- 
batic (Born-Oppenheimer) approximation, which allows 
us to define a potential energy surface for the nuclear 
coordinates. As in classical molecular dynamics sim- 



ulations, an important issue of the PIMD method is 
the proper description of interatomic interactions, which 
should be as realistic as possible. To overcome the 
limitations of effective classical potentials to reproduce 
the many-body energy surface, one has to resort to 
self-consistent quantum-mechanical methods. Neverthe- 
less, true density functional or Hartree-Fock based self- 
consistent potentials require computer resources that 
would restrict enormously the size of our simulation cell. 
We have found a reasonable compromise by deriving the 
Born-Oppenheimer surface for the nuclear dynamics from 
an efficient tight-binding effective Hamiltonian, based on 
density functional (DF) calculationsj^ The capability of 
tight-binding methods to simulate different properties of 
solids and molecules has been reviewed by Goringe et 
aZ.'^^ We have checked the abihty of this DF-TB potential 
to predict frequencies of C-H vibrations. In particular, 
for a methane molecule it predicts in a harmonic approxi- 
mation frequencies of 3100 and 3242 cm~^ for C-H modes 
with symmetry Ai and T2, respectively, to be compared 
with experimental frequencies of 2917 and 3019 cm~^iM 
Taking into account the typical anharmonic shift asso- 
ciated to these modes (towards lower frequencies), the 
agreement is satisfactory. A detailed analysis of vibra- 
tional frequencies of hydrocarbon molecules derived with 
the present DF-TB potential, including anharmonicities, 
can be found elsewhere. ■^^•'^^ 

Simulations were carried out on a 2 x 2 x 2 super- 
cell of the diamond face-centered cubic cell with periodic 
boundary conditions, containing 64 C atoms and one im- 
purity. For comparison, we also carried out simulations 
of diamond without impurities, using the same supercell 
size. In our calculations, four symmetry independent k 
points in the Brillouin zone of the simulation supercell 
were employed, distributed by following the Monkhorst- 
Pack generation schemei^i We have checked the conver- 
gence of the internal energy for some selected atom con- 
figurations, by considering up to 32 k points. In partic- 
ular, the internal energy of hydrogen on the tetrahedral 
T site and on the BC configuration changes by about 1 
meV, respect to our calculation employing 4 k points. 
Sampling of the configuration space has been carried out 
at temperatures between 100 and 800 K. The simulation- 
cell parameter employed in our calculations was taken 
from experimental data^ and thus changed from 7.1336 
A at 100 K to 7.1424 A at 800 K. For a given tempera- 
ture, a typical run consisted of 10^ MD steps for system 
equilibration, followed by 10^ steps for the calculation of 
ensemble average properties. 

To have a nearly constant precision in the path inte- 
gral results at different temperatures, we have taken a 
Trotter number that scales as the inverse temperature. 
At 300 K we took L = 20 for H and D, and L = 70 
for Mu. For comparison with the results of our PIMD 
simulations, we have carried out some classical MD sim- 
ulations with the same interatomic interaction (setting L 
= 1). Note that for equilibrium properties in the canoni- 
cal ensemble (i.e., mean energy, spatial distribution), the 
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classical limit is equivalent in our context to the infinite- 
mass limit (m — *■ oo). The quantum simulations were 
performed using a staging transformation for the bead co- 
ordinates. Chains of four Nose-Hoover thermostats with 
mass Q — f3h'^/5L were coupled to each degree of free- 
dom to generate the canonical ensemblei^ To integrate 
the equations of motion, we have used the reversible ref- 
erence system propagator algorithm (RESPA), which al- 
lows one to define different time steps for the integra- 
tion of the fast and slow degrees of freedomrS*^ The time 
step At associated to the calculation of DF-TB forces 
was taken in the range between 0.2 and 0.5 fs, which 
was found to be appropriate for the interactions, atomic 
masses, and temperatures considered here. For the evo- 
lution of the fast dynamical variables, that include the 
thermostats and harmonic bead interactions, we used a 
smaller time step 5t = At /A. We note that for muonium 
at 200 K (the lowest temperature considered for this im- 
purity), a simulation run of 10'^ MD steps requires the 
calculation of forces and energy with the tight-binding 
code for about 10'' configurations, which needs the use of 
large parallel computers. 



B. Calculation of anharmonic vibrational 
frequencies 

Important characteristics of impurities in solids are 
their vibrational frequencies, which are known to be de- 
pendent on the particular site occupied by a given impu- 
rity in the crystal and on its interactions with the nearby 
hosts atoms. In our context, the question arises whether 
the oscillator frequencies associated to an impurity lo- 
cated at a given site can be extracted by assuming the 
host C atoms fixed in the relaxed geometry. This is in 
fact a method frequently employed to calculate vibra- 
tional frequencies of impurities in solids. On the other 
side, when all C atoms are allowed to relax by follow- 
ing the impurity motion, the potential energy surface is 
much flatter (smaller energy changes) than when the host 
atoms are fixed. To obtain an approach for the actual vi- 
brational frequencies of the impurities, one can calculate 
the eigenvalues of the dynamical matrix of the whole sim- 
ulation cell, and obtain the frequencies in the harmonic 
approximation (HA). However, for light impurities and 
atomic configurations displaying large strains (important 
relaxations), the anharmonicity can be appreciable, mak- 
ing the harmonic frequencies meaningless. 

To calculate these anharmonic frequencies we will em- 
ploy a method based on the linear response (LR) of the 
system to vanishingly small forces applied on the atomic 
nuclei. With this purpose, we consider a LR function, the 
static isothermal susceptibility i that is readily derived 
from PIMD simulations of the equilibrium solid, without 
dealing explicitly with any external forces in the simula- 
tion. This approach represents a significant improvement 
as compared to a standard harmonic approximation^^ A 
sketch of the method is given in the following. 



Let us call {xip} the set of 3NL Cartesian coordinates 
of the beads forming the ring polymers in the simula- 
tion cell {i = 1, . . . , 3N;p = 1, . . . , L). We consider the 
set {Xi} of centroid coordinates, with Xi defined as the 
mean value of coordinate i over the corresponding poly- 
mer: 

1 ^ 

Xi = - "^x^p . (1) 

Then, the linear response of the quantum system to small 
external forces on the atomic nuclei is given by the sus- 
ceptibility tensor x"^, that can be defined in terms of 
centroid coordinates aa^ 

xlj = (i^fm~m~ fi^j , (2) 

where /? = {kBT)~^, rrii is the mass of the nucleus asso- 
ciated to coordinate i, fiij = (Xa,) - {Xi){Xj) is the 
covariance of the centroid coordinates Xi and Xj, and 
(...) indicates an ensemble average along a MD run. 

The tensor allows us to derive a LR approximation 
to the low-lying excitation energies of the vibrational sys- 
tem, that is applicable even to highly anharmonic situ- 
ations. The LR approximation for the vibrational fre- 
quencies reads 



where A„ {n ~ 1,...,3A'^) are eigenvalues of x"^, and 
the LR approximation to the low-lying excitation energy 
of vibrational mode n is given by huJn,LR- More details 
on the method and illustrations of its ability for predict- 
ing vibrational frequencies of solids and molecules can 
be found elsewhere i^^i^Si^ii^ In connection with the vi- 
brational modes that actually appear in our calculations, 
we note that the application of periodic boundary condi- 
tions is physically equivalent to the only consideration of 
lattice vibrations at the center {k = 0) of the Brillouin 
zone of the simulation celli^ 



III. RESULTS 
A. Energy 

We first discuss the stable sites for the hydrogenic im- 
purities in diamond, as derived from classical calculations 
at T = 0, i.e., point nuclei without spatial delocalization. 
In this respect, Estle et ali noticed the importance of 
lattice relaxation for finding the most stable site for the 
impurities, and obtained that interstitial bond-centered 
hydrogen or muonium is stable as a result of unusually 
large lattice relaxation. This is in fact the case for the 
DF-TB model employed here. After relaxation of the 
host atoms, the energy surface has minima only at two 
non-equivalent positions. The lowest-energy position is 
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the BC site, as in several earlier investigationsii^ In the 
relaxed configuration, we find a distance C-H of 1.17 A, 
which means a backward relaxation of the adjacent C 
atoms of 0.40 A. This indicates that the formation of the 
C-H-C bridge requires a large dilation of the C-C bond 
(a 52% of the original one). These values agree with those 
found earlier from tight-binding calculations of hydrogen 
in diamond. 

Another local minimum of the energy surface is found 
for the impurity at the tetrahedral T site, located 1.44 eV 
above the absolute minimum. In this configuration, the 
relaxation of the nearest lattice atoms is much smaller 
than for the bond-center configuration (0.08 A). There 
appears in the literature a rather large dispersion of val- 
ues for the relative energy of H on the T site. These 
values range typically from 0.5 to 2.7 eV»i^ Then, our re- 
sult lies in the intermediate region of these earlier values, 
and in particular it is close to that reported by Kauko- 
nen et al*^ from a DF-TB calculation (1.6 eV). For the 
energy barrier from a BC to a T site we find 1.7 eV, 
a value between those obtained from local-density func- 
tional theorjiiS (1.6 eV) and earlier DF-TB calculations'*'^ 
(2.0 ±0.1 eV). For the migration barrier from T to BC 
sites we have found 0.3 eV, to be compared with 0.4±0.1 
eV obtained in Ref. ^3- We note that the so-called anti- 
bonding site, between an atom and a T site along a (111) 
direction, has turned out to be from our TB calculations 
a saddle point on the energy surface, contrary to Ref. 0, 
where a local minimum was reported. 

We now turn to our simulations at finite temperatures. 
The internal energy, E{V,T), at volume V and tempera- 
ture T can be written as: 



E{V,T) = En,in{V) + E,{V,T) , 



(4) 



where -Emin(^) is the potential energy for the classical 
solid at T = 0, and Ev{V,T) is the vibrational energy. 
At finite temperatures, V changes with T due to thermal 
expansion. The vibrational energy, E^{V,T), depends 
explicitly on both V and T, and can be obtained by sub- 
tracting the energy £'min(^) from the internal energy. In 
this way, path-integral molecular dynamics simulations 
allow us to obtain E^{V, T) for a given volume and tem- 
perature. 

Shown in Fig.^ is the temperature dependence of 
the vibrational energy E^ per simulation cell for hydro- 
genic impurities at the BC site. Symbols indicate re- 
sults of PIMD simulations: deuterium (circles), hydrogen 
(squares), and muonium (triangles). For comparison we 
also show iJv for pure diamond (triangles). At 300 K, 
the vibrational energy of diamond amounts to 12.4 eV 
per simulation cell, i.e., 0.19 eV/atom. As expected, E^ 
increases as temperature is raised, and eventually con- 
verges to the classical limit £'^' = SNksT at high T. 
When we consider the hydrogenic impurities at the BC 
site, the vibrational energy increases as compared to pure 
diamond. This increase is larger the lower the impurity 
mass, as a consequence of zero-point motion. 
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FIG. 1: Temperature dependence of the vibrational energy of 
the 2x2x2 diamond supercell with one impurity, as derived 
from PIMD simulations. Results are shown for deuterium 
(circles), hydrogen (squares), and muonium (triangles). For 
comparison, we also present results for a pure diamond su- 
percell (diamonds). Dashed lines are guides to the eye. Error 
bars are less than the symbol size. 



At this point, an interesting characteristic of the dif- 
ferent hydrogenic defects is their associated vibrational 
energy. At a given temperature, this energy is defined as 
the difference 1\E^ = E^{ioAC -I- Imp) - £;v(64C), where 
'Imp' stands for H, D, or muonium. Note that Ai?v so 
defined is not just the vibrational energy of a given im- 
purity in the crystal, but it also includes changes in the 
vibrations of the host atoms. In Fig.[21 the energy Ai5v 
is plotted vs temperature for the considered impurities 
at a BC site. Within the numerical precision of our MD 
simulations, Ai^v turns out to be constant for the three 
hydrogenic impurities. This means that at temperatures 
lower than 800 K we find for these defects, Ai?v values ba- 
sically the same as their corresponding zero-temperature 
limit. A slight increase of this energy should be expected 
in the temperature region shown in Fig.|21 especially for H 
and D (with vibrational frequencies lower than for Mu), 
due to thermal excitation of higher vibrational modes. 
However, it seems that this increase is compensated for 
by a decrease due to lattice expansion, with its associ- 
ated softening of the vibrational modes. In fact, we have 
carried out a PIMD simulation of H at a BC site with the 
lattice parameter corresponding to 200 K (a = 3.5668 A), 
and found AiJv = 0.40 eV, about 50 meV higher than 
the value obtained using the actual lattice parameter at 
800 K (a = 3.5712 A). For comparison with our results 
of PIMD, we also show in Fig. |21 results for A_Ev derived 
from classical MD simulations (diamonds). These data 
points follow closely the trend expected for a classical 
three-dimensional harmonic oscillator: Ai?v — SkBT. 
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TABLE I: Vibrational energy A_Ev of liydrogenic defects and position of the defect level, Ei, at 300 K for impurities at BC 
and T sites, as derived from PIMD simulations. Energy is given in eV, and the zero of energy for the electronic levels is taken 
at the valence-band top. Also given are the ratios between vibrational energies for different impurities. Error bars in AE^ are 
± 6 meV for H and D, and ± 10 meV for Mu. The statistical error for Ei is ±2 meV . 
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FIG. 2: Vibrational energy of the hydrogenic defects as a 
function of temperature. Results are shown for deuterium 
(circles), hydrogen (squares), and muonium (triangles). For 
comparison, we also present the vibrational energy derived 
from classical MD simulations. Error bars are on the order of 
the symbol size. Dotted lines are guides to the eye. 



In a one-particle harmonic approximation, the low- 
temperature vibrational energies corresponding to deu- 
terium, hydrogen, and muonium scale as 0.71 : 1 : 2.97, 
i.e., proportional to rnT^I"^ (rn, impurity mass). In Table 
I we give the vibrational energies of the hydrogenic de- 
fects and ratios between them. At 300 K, Ai?v values de- 
rived from our PIMD simulations scale as 0.75 : 1 : 2.86, 
with ratios somewhat different from those for a har- 
monic approximation. In fact, we have estimated error 



bars of ±0.02 and ±0.05 for the ratios IAE° j AE^ and 
AE^^ / AE^ , respectively. This is an indication of the 
anharmonicity present in the defect vibrations, which will 
be discussed below. 

By using the same procedure, we have calculated the 
vibrational energy AE^ for impurities at the T site. The 
results at 300 K are given in Table I. The resulting en- 
ergies are larger than the corresponding ones for the BC 
configuration. At first sight, this result may seem surpris- 
ing, taking into account the constrained geometry of the 
BC configuration, which should give a larger frequency 
for the stretching vibration along the C-H-C axis. How- 
ever, vibrations transverse to this axis are expected to 
have much lower frequency, making reasonable a larger 
Ai?v for the T defect with threefold degenerate vibra- 
tional modes (see below). For the energy ratios at the T 
configuration we find 0.72 : 1 : 3.01, which are compatible 
within error bars with a harmonic approach. 

A detailed analysis of the temperature dependence of 
properties of the impurities at the T site is not feasible 
due to the metastability of these defects. The metastabil- 
ity of hydrogen at the T site is illustrated in Fig.O where 
we show the energy of the simulation cell along a PIMD 
run at 600 K. This simulation run started with H at a T 
site, and after about 15000 steps we observed the passage 
of the impurity to a BC site. This is reflected in the figure 
by an energy jump of about 1.6 eV, which corresponds 
to the difference between energy minima plus the differ- 
ence between vibrational energies. Similar results were 
obtained in other simulations at the same temperature, 
with hydrogen typically jumping from T to BC positions 
after a number of MD steps in the range from 10** to 
10^. However, at 300 K hydrogen remained around the 
metastable T site during several simulation runs of 10^ 
steps, and the same happened for deuterium and muo- 
nium. We remember that the time scale of a simula- 
tion run (i.e.. At times the number of MD steps) does 
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FIG. 3: Internal energy of the simulation cell along a PIMD 
trajectory at 600 K, starting with hydrogen on a T site. One 
observes a drop in the energy after about 15000 steps, cor- 
responding to a passage of hydrogen from the T-site region 
to the BC well. The zero of energy is taken for a classical 
system with the impurity on a BC site at T = (absolute 
energy minimum). 



not represent a real physical time, since the dynamics of 
the beads in the ring polymers does not correspond to 
the actual dynamics of the quantum particles. The fic- 
titious bead dynamics is merely used as an efficient way 
to calculate equilibrium properties of the system in the 
canonical ensemble. Thus, the plot in Fig. 13 illustrates 
the ability of PIMD to explore the configuration space, 
and to drive the system to the most stable region, once 
the temperature is high enough. For comparison, we note 
that the metastability of muonium at the T site has been 
observed by muon spin rotation44 In fact, 'isotropic muo- 
nium' (at T) converts to 'anisotropic muonium' (at BC) 
in the range 500-700 K. 



B. Vibrational frequencies 

As indicated above, the DF-TB potential gives a 
good description of vibrational modes in hydrocarbon 
molecules. To check the reliability of this potential to 
give vibrational frequencies of hydrogenic impurities in 
diamond, and to assess their anharmonicity, we have 
carried out calculations for a simple atomic cluster (see 
Fig.Q. In this cluster, the impurity occupies a position 
similar to a BC site in the diamond lattice, midway be- 
tween two carbon atoms with a distance d(C,H) = 1.17 A, 
as that found in diamond with the DF-TB potential. In 
this atomic cluster, we have calculated total energies for 
hydrogen displacements along the C-C axis, with both 
our TB Hamiltonian and the B3LYP density-functional 




B3LYP, COha= 



FIG. 4: Schematic representation of the atomic cluster em- 
ployed to compare total energy changes obtained with our 
tight-binding approach and with the B3LYP model, ujha in- 
dicates the vibrational frequency for hydrogen along the C-C 
direction, obtained in a single-particle harmonic approxima- 
tion in either case. 



theory (with the cc-pDVZ basis set as implemented in the 
Gaussian 98 package). Results obtained by both meth- 
ods are shown in Fig. [51 where the zero of energy corre- 
sponds to the BC position in either case. The continu- 
ous and dashed lines correspond to DF-TB and B3LYP 
models, respectively. The curvature of the solid line at 
the minimum is smaller than that of the dashed line, and 
thus the harmonic approach gives for the DF-TB method 
a vibrational frequency (1991 cm~^ for H) lower, but of 
the same order, than that of the B3LYP potential (2293 
cm~^ for H). We emphasize that these energy curves do 
not correspond to hydrogenic impurities in bulk diamond, 
but they help to analyze the reliability of the DF-TB 
potential and to illustrate qualitatively the vibrational 
potential of these impurities in the solid. 

In Fig. Owe also present the two lowest energy levels 
for hydrogen, obtained by numerically solving the one- 
dimensional Schrodinger equation with each interaction 
model. Now the ground state for the B3LYP potential 
is higher than that of the DF-TB model. The first ex- 
cited states yielded by both procedures are very close 
to one another, and in fact are indistinguishable on the 
scale of Fig.O We then find that the excitation energy 
h{uji —oJo) corresponds to 2633 and 2577 cm~^ for DF-TB 
and B3LYP, respectively. In both cases we find excita- 
tion energies larger than in the harmonic approximation, 
i.e., anharmonicity hardens the stretching vibration, as 
expected for an impurity in a highly-confined geometry. 
The excitation energies obtained by both methods are 
similar, with a relative difference between them smaller 
than 3%, which gives us confidence in the reliability of 
the DF-TB method to analyze vibrational properties of 
hydrogenic point defects in diamond. 

We now go back to the study of hydrogenic impuri- 
ties in bulk diamond. Before analyzing the anharmonic 
modes given by the linear response procedure described 
in Sec. II. B, we will present the vibrational frequencies 
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FIG. 5: Potential energy curves for impurity displacements 
along the direction parallel to the C-C axis, with fixed carbon 
atoms, in the cluster shown in Fig.^] The minimum potential 
energy corresponds to the impurity at the bond-center site. 
Horizontal lines represent energy levels for hydrogen, obtained 
by numerically solving the one-dimensional Schrodinger equa- 
tion with those potentials. Results are shown for the tight- 
binding potential employed here (solid lines) and for the 
B3LYP model (dashed lines). The given frequencies corre- 
spond to first excitation energies in both cases. 



of the impurities derived from a harmonic approxima- 
tion. This will allow us to assess the importance of an- 
harmonicity for the impurity vibrations. By calculating 
the dynamical matrix, we find the vibrational density 
of states (VDOS) corresponding to the whole simulation 
cell (64 C atoms plus one impurity). We then use the 
amplitude of the displacement of the impurity in each 
vibrational mode, to obtain its partial VDOS. This is 
shown in Fig.|Hlfor (a) deuterium, (b) hydrogen, and (c) 
muonium at the BC site. The results can be most clearly 
understood for the case of muonium, where one observes 
basically two peaks at 2221 and 4775 cm~^, well sepa- 
rated from the region of crystal vibrations (up to about 
1500 cm"-'^). Due to the low mass of Mu, these vibrations 
are almost totally decoupled from the host-atom vibra- 
tions, and are associated to motion perpendicular (bend- 
ing, twofold degenerate) and parallel (stretching) to the 
C-C axis, respectively. This is in fact the situation ex- 
pected for an impurity at a BC site, assuming fixed host 
atoms. The partial VDOS is, however, more complicated 
for H and D, which participate in modes that are not 
totally localized at the impurity atom or at the impurity 
plus its nearest C neighbors. For H, we find only one 
mode with frequency larger than the lattice modes, cor- 
responding to the stretching vibration along the C-H-C 
bond (w||,ffA = 1738 cm~^). We obtain two other peaks 
at 448 and 719 cm~^ with large participation of hydro- 
gen motion. They are assigned to transverse vibrations, 
perpendicular to the C-H-C axis, and include impor- 
tant participation of nearby C atoms. For deuterium, 
the most prominent feature is a large peak at 367 cm~^, 
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FIG. 6: Vibrational density of states for impurities at a BC 
site, as derived from the dynamical matrix of the simulation 
cell in a harmonic approximation. Each panel corresponds 
to an impurity: (a) deuterium, (b) hydrogen, (c) muonium. 
Labels indicate wave numbers in cm~^. For the sake of clarity, 
discrete modes have been broadened to Gaussian profiles with 
a standard deviation of 10 cm~^. Note the different horizontal 
scale in panel (c). 



corresponding to bending vibrations. In this case, the 
stretching vibration is highly coupled to crystal modes, 
and does not appear as a prominent peak in the VDOS. 

Next, we apply the linear-response procedure to obtain 
the (anharmonic) vibrational frequencies of the impuri- 
ties. Again, we obtain the VDOS of the whole simulation 
cell and project from it the contribution of the considered 
impurity. This partial VDOS is shown in Fig.[7| for the 
hydrogenic impurities studied here at 300 K. As before, 
we first comment on the results for muonium in panel (c). 
We observe two distinct peaks, as for the HA, but now 
at much higher frequencies. This resembles the effect 
of anharmonicity observed for the atomic cluster con- 
sidered above. For muonium at a BC site in diamond, 
both vibrations parallel and perpendicular to the C-Mu- 
C axis harden strongly due to anharmonicity. Something 
similar happens for hydrogen and deuterium, as noticed 
when one compares the corresponding panels in Figs. El 
and For hydrogen, an interesting feature is the ap- 
pearance of a second peak at 1559 cm"'^, above but close 



to the largest frequency in bulk diamond, corresponding 
to a bending mode [Fig.[7{b)]. For deuterium, we ob- 
tain a well-resolved stretching mode at 1837 cm~^, far in 
frequency from the lattice vibrations. For the frequen- 
cies derived from the LR method, we estimate error bars 
of 20 cm^^, due to the statistical noise present in the 
susceptibility tensor derived from the PIMD simula- 
tions. For the different impurities, the highest-frequency 
mode 1^11, L_R scales as 0.72 : 1 : 3.13, to be compared 
with the ratio 0.71 : 1 : 2.97 expected for one-particle 
harmonic oscillators. The stretching frequency derived 
for H from our LR calculations {oj\\^lr = 2544 cm~^) 
is lower than that obtained for this impurity by Goss et 
alm^ii^ from DF calculations in the harmonic approxima- 
tion (2919 cm-i). This difference can be explained, at 
least in part, by the difference in the C-H distance ob- 
tained by both methods: 1.17 A for our TB potential vs 
1.13 A for the DF calculations^^ (the larger the distance, 
the lower is expected to be this frequency) . We also note 
that DF calculations predict a frequency of 2456 cm^^ for 
positively charged hydrogen in a buckled bond-centered 
configurationpi^ 

We have repeated the same LR procedure to calcu- 
late the vibrational frequencies of hydrogenic impurities 
at the T site. In this case, we find for each impurity a 
threefold degenerate mode, in agreement with the tetra- 
hedral point symmetry of this site. Coupling with lattice 
vibrational modes is negligible. We find the frequencies 
2088, 2818, and 7284 cm"! for D, H, and Mu, respec- 
tively, which translates to a frequency ratio 0.74:1:2.58. 
For impurities at the T site, anharmonicity causes a soft- 
ening of the modes. For example, for hydrogen we find in 
the HA a frequency of 3099 cm^^, more than 200 cm~^ 
larger than that found with the LR method. This soften- 
ing is still more important for muonium, due to its larger 
zero-point vibration, and is the origin of the small ratio 
^\\lr/^\\ LR ~ ^-^^ compared with a ratio of 2.97 

expected for harmonic oscillators). We note that anhar- 
monicity affects the mode frequencies in opposite ways 
at the BC and T sites. While the vibrational modes as- 
sociated to the impurities at the BC site suffer a strong 
shift towards higher frequencies, for the impurities at the 
tetrahedral site the shift is towards lower ones. This no 
doubt is related to the different geometries in both con- 
figurations: at the BC site the impurity is in a much more 
confined environment, with a strongly directional char- 
acter (along the bond), while at the T site the impurity 
is less constrained. 



C. Defect levels in the electronic gap 

The DF-TB method employed here allows us to study 
the influence of zero-point motion on the one-electron 
levels. For a classical hydrogen-like impurity at a BC 
site, we find a defect level located in the electronic gap 
2.61 eV above the top of the valence band. This result 
is close to the donor level found by Goss et al. from 
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FIG. 7: Vibrational density of states for impurities at a BC 
site and T — 300 K, as derived from a linear-response analysis 
of PIMD results. Each panel corresponds to an impurity: (a) 
deuterium, (b) hydrogen, (c) muonium. Labels indicate wave 
numbers in cm~^. For clarity, discrete modes are presented 
with a Gaussian profile with a standard deviation of 10 cm~^. 
Note the different horizontal scale in panel (c). 



local-density- functional calculations.'^^ Also, Saada et al. 
found a defect level at ~ 0.5 eV above the middle of 
the energy gap for the minimum-energy configuration, 
which in their case corresponded to the impurity on off- 
BC sitesii 

Defect levels associated to the presence of the impurity 
in the crystal are expected to change as a consequence of 
electron-phonon interaction. In fact, effects of this kind 
of interactions on the electronic structure have been ex- 
perimentally investigated by measuring the temperature 
dependence of the optical spectra of solidsi^ These ef- 
fects appear even at low T, due to zero-point motion, and 
can be large in some cases. For example, in diamond the 
direct electronic gap (at the F point) is reduced by about 
0.6 eV with respect to that calculated when zero-point 
vibrations are neglected4& Such studies have been com- 
plemented in some cases by observing the dependence of 
the spectra on isotopic mass, whenever different stable 
isotopes are available424& 

Similar effects of quantum atomic vibrations should 
also appear in the defect levels induced by the hydrogenic 



9 



0.05 



-1 1 1 1 1 1 r- 



classical 



-0.05 



^ -0.10 



-0.15 



D 



H 



Mu 



For hydrogenic atoms at site T, the level goes up as 
the impurity mass is lowered, contrary to the case of the 
bond-center configuration. Also, the level for the T con- 
figuration is more sensitive to the impurity mass, since it 
changes by 0.47 eV from a classical impurity to muonium, 
vs a decrease of 0.10 eV for the BC configuration. 

We note that, while the absolute position of the de- 
fect level can depend on the model employed, its rela- 
tive change induced by temperature and isotopic mass 
is expected to be quite reliable, because in our analysis 
both electronic and nuclear dynamics are treated at the 
level of quantum mechanics, which permits a more realis- 
tic description of the electron-phonon interaction in this 
system. 
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FIG. 8: Temperature dependence of the defect-level energy 
for deuterium (circles), hydrogen (squares), and muonium 
(triangles) at a bond center. For comparison, we also present 
results for a classical simulation (diamonds). The zero of en- 
ergy corresponds to the level of a classical impurity at T = 0. 
Error bars are on the order of the symbol size. Dashed lines 
are guides to the eye. 

impurities studied here. Shown in Fig.|Hlis the tempera- 
ture dependence of the defect level energy for H, D, and 
Mu at a BC site in diamond, as derived from our PIMD 
simulations. For comparison, we also present the energy 
level Ej derived from classical MD simulations. In this 
Figure, the zero of energy corresponds to the classical 
limit at T = 0. The most remarkable fact that one ob- 
serves in this plot is an appreciable shift of the energy 
level as the impurity mass changes. In particular, for 
the impurity at the BC site, the level goes down as the 
mass is lowered. At 300 K we find a decrease in Ej of 13 
meV when passing from D to H, and a further decrease 
of 70 meV when comparing muonium with deuterium. 
We note also that the level in the case of deuterium is 
20 meV lower than that corresponding to a classical im- 
purity (limit of infinite mass) at 300 K. In all cases, the 
statistical error due to finite sampling of the canonical 
ensemble is of about 2 meV. 

For a given impurity, the defect level goes down as 
the temperature is raised. In fact, for H and D we find 
a decrease in the level energy of about 30 meV when 
T is increased from 100 to 800 K. For muonium, this 
temperature-induced change amounts to about 20 meV. 
For the classical impurity this level shift is more impor- 
tant, resulting to be 63 meV. 

The position of the defect level for the different impu- 
rities at 300 K is given in Table II, where we also present 
the energy of the level for impurities at the tetrahedral 
T site. In this case, the classical result at 300 K is found 
to be 0.58 eV above the defect level for the impurity at 
a BC site. 



IV. DISCUSSION 

In Sec. Ill we have presented results of our PIMD sim- 
ulations for H, D, and Mu in diamond. The main advan- 
tage of this kind of simulations is the possibility of calcu- 
lating defect energies at finite temperatures, with the in- 
clusion of a full quantization of host-atom motions, which 
are not easy to be accounted for in fixed-lattice calcula- 
tions. Isotope effects can be readily explored, since the 
impurity mass appears as an input parameter in the cal- 
culations. This includes the consideration of zero-point 
motion, which together with anharmonicity causes ap- 
preciable non-trivial effects. Our results indicate that 
hydrogenic impurities at the bond-center site cannot be 
accurately described as a particle moving in a harmonic 
potential. Even if anharmonicities of the interatomic po- 
tential are taken into account, a single-particle approx- 
imation does not give reliable results for the impurity 
complex at finite temperatures. It is then necessary to 
treat the defect as a many-body problem with anhar- 
monic interactions. 

The most prominent feature in the VDOS of the hy- 
drogenic impurities at a BC site is a peak corresponding 
to the stretching vibration, with a frequency higher 
than the lattice vibrations. This stretching mode is 
infrared-active, but to our knowledge has not been de- 
tected yet. On the theoretical side, there are not many 
calculations in the literature predicting the frequency of 
this mode. Goss et alm^ employed an HA from local- 
density-functional calculations and found for hydrogen 
at a BC site (in the neutral charge state) a stretching 
frequency uj\\ = 2919 cm^-'^, somewhat higher than our 
result using LR calculations Li?, = 2544 cm"-'^). 

From PIMD simulations we have calculated in Sec. 
III. A the defect energy AE^ for the different hydrogenic 
impurities. Also, by using the LR procedure we have ob- 
tained an approach to the vibrational frequencies of the 
impurities, which is more realistic than a pure harmonic 
approximation. These frequencies can now be used to 
estimate a zero-point energy for the impurities in dia- 
mond. In the case of muonium, the modes are almost 
totally localized on the impurity, and the zero-point vi- 
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brational energy of the defect is well approximated by 
^K,LR = + ^||,Lfl/2). This gives AE° ,^i^ = 

0.96 eV, close to the ground state energy of the defect 
presented in Fig.^ (0.99 ± 0.01 eV). Using the frequen- 
cies uj± HA and cjii HA derived from the HA, one finds 
HA = 0.57 eV, much lower than that obtained di- 
rectly from PIMD simulations. Note that we have used 
here an expression for Eq which is only rigorously valid 
for harmonic oscillators, but it nevertheless gives a good 
approximation to the zero-point energy of the defect 
when one introduces the frequencies derived from the 
LR method. In this sense, these LR frequencies may be 
considered as renormalized harmonic frequencies, which 
incorporate in a nonperturbative way anharmonicities of 
the interatomic interactions, as seen by the atoms in their 
quantum motion. 

The zero-point energy of the defect requires a com- 
ment in the sense that, in the case of muonium, it could 
help to cross the adiabatic barrier for impurity diffusion. 
Nevertheless, we find that muonium is confined in the 
BC region, especially at low temperatures. Impurity mi- 
gration to other regions of the crystal requires important 
relaxations of the neighboring C atoms, which are not 
probable at low temperatures. This picture is similar to 
that described in the literature as "opening of a door" ^ 
which favors impurity diffusion. 

An analysis of hydrogen diffusion in diamond is out of 
the scope of this paper. As noted above, actual diffusion 
coefficients are not directly accessible with the kind of 
MD simulations employed here, since the time scale in 
the simulations is not readily connected to the real one. 
In connection with this, PIMD simulations can be applied 
to study quantum diffusion (including tunneling) of hy- 
drogen in pure and doped diamond, by calculating free- 
energy barriers in a way similar to that employed earlier 
to study H diffusion in metals^^ and semiconductors!^ 
Even though the rate of H tunneling in diamond is not 
expected to be high at low T, due to the large lattice 
relaxation associated to the impurity on a BC site, ther- 
mally activated tunneling may be possible, as observed 
for hydrogen in siliconi^ This point will require further 
research. 

Theoretical techniques to calculate the electronic band 
structure of solids have been improving their precision for 
many years. For various purposes, the accuracy currently 
achieved by these methods is excellent, when comparing 
their predictions with experimental data. However, zero- 
point motion is a significant factor limiting the accuracy 
of state-of-the-art techniques to predict energy bands. 



The same happens for defect levels caused by impurities 
in solids, since their energy may change appreciably as 
the impurity mass is varied. This effect has been illus- 
trated here by shifts in the energy of the levels due to 
hydrogenic impurities in BC or T sites in diamond. For 
muonium at the T site, we have found a level shift of 0.47 
eV (upwards), on the order of the renormalization of the 
diamond gap due to zero-point motion (^ 0.6 eV). For 
the BC site, the corresponding shift for Mu was 0.10 eV 
and for H we found 0.03 eV (both downwards). Thus the 
magnitude and direction of these shifts depend markedly 
on the mass and position of the impurity in the solid. 
This behavior can be rationalized by considering the re- 
sults of perturbation theory for the electron-phonon in- 
teraction in the limiting case of a harmonic oscillator 
The shift of an electronic level is expected to be propor- 
tional to the square of the mode amplitude, i.e., propor- 
tional to mT^/"^ in the low-temperature limit and to T at 
high temperatures. The extrapolation to zero tempera- 
ture of the defect levels for D, H, and Mu at a BC site (see 
Fig.© gives shifts scaling as 0.72 : 1 : 2.24, to be com- 
pared to the ratio 0.71 : 1 : 2.97 expected for a harmonic 
impurity. The energy shift for Mu deviates appreciably 
from the harmonic expectation, as a consequence of the 
larger anharmonicity of this impurity center. 

In summary, the PIMD method has turned out to be 
well-suited to study finite-temperature equilibrium prop- 
erties of light impurities in diamond. This has allowed 
us to notice the importance of anharmonicity in order 
to give a realistic description of the properties of these 
point defects. This anharmonicity shows up in the vibra- 
tional modes of the impurities, causing important shifts 
respect to the harmonic expectancy. Also, zero-point mo- 
tion introduces an appreciable shift in the defect levels, 
which depends on the impurity mass. This type of analy- 
sis offers a promising way for studying other challenging 
effects of light impurities in diamond, such as quantum 
diffusion. 
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